Extensions to the Proximal Distance Method of Constrained Optimization

The current paper studies the problem of minimizing a loss f(x) subject to constraints of the form Dx ∈ S, where S is a closed set, convex or not, and D is a matrix that fuses parameters. Fusion constraints can capture smoothness, sparsity, or more general constraint patterns. To tackle this generic class of problems, we combine the Beltrami-Courant penalty method of optimization with the proximal distance principle. The latter is driven by minimization of penalized objectives f(x)+ρ2dist(Dx,S)2 involving large tuning constants ρ and the squared Euclidean distance of Dx from S. The next iterate xn+1 of the corresponding proximal distance algorithm is constructed from the current iterate xn by minimizing the majorizing surrogate function f(x)+ρ2‖Dx−𝒫S(Dxn)‖2. For fixed ρ and a subanalytic loss f(x) and a subanalytic constraint set S, we prove convergence to a stationary point. Under stronger assumptions, we provide convergence rates and demonstrate linear local convergence. We also construct a steepest descent (SD) variant to avoid costly linear system solves. To benchmark our algorithms, we compare their results to those delivered by the alternating direction method of multipliers (ADMM). Our extensive numerical tests include problems on metric projection, convex regression, convex clustering, total variation image denoising, and projection of a matrix to a good condition number. These experiments demonstrate the superior speed and acceptable accuracy of our steepest variant on high-dimensional problems. Julia code to replicate all of our experiments can be found at https://github.com/alanderos91/ProximalDistanceAlgorithms.jl


Introduction
The generic problem of minimizing a continuous function f (x) over a closed set S of R p can be attacked by a combination of the penalty method and distance majorization.The classical penalty method seeks the solution of a penalized version h ρ (x) = f (x) + ρq(x) of f (x), where the penalty q(x) is nonnegative and 0 precisely when x ∈ S. If one follows the solution vector x ρ as ρ tends to ∞, then in the limit one recovers the constrained solution (Beltrami, 1970;Courant, 1943).The function is one of the most fruitful penalties in this setting.Our previous research for solving this penalized minimization problem has focused on an MM (majorization-minimization) algorithm based on distance majorization (Chi et al., 2014;Keys et al., 2019).In distance majorization one constructs the surrogate function using the Euclidean projection P(x n ) of the current iterate x n onto S. The minimum of the surrogate occurs at the proximal point According to the MM principle, this choice of x n+1 decreases g ρ (x | x n ) and hence the objective h ρ (x) as well.As we note in our previous JMLR paper ( Keys et al., 2019), the update (1) reduces to the classical proximal gradient method when S is convex (Parikh, 2014).
We have named this iterative scheme the proximal distance algorithm (Keys et al., 2019;Lange, 2016).It enjoys several virtues.First, it allows one to exploit the extensive body of results on proximal maps and projections.Second, it does not demand that the constraint set S be convex.If S is merely closed, then the map P(x) may be multivalued, and one must choose a representative element from the projection P(x n ).Third, the algorithm does not require the objective function f (x) to be differentiable.Fourth, the algorithm dispenses with the chore of choosing a step length.Fifth, if sparsity is desirable, then the sparsity level can be directly specified rather than implicitly determined by the tuning parameter of the lasso or other penalty.
Traditional penalty methods have been criticized for their numerical instability.This hazard is mitigated in the proximal distance algorithm by its reliance on proximal maps, which usually are highly accurate.The major defect of the proximal distance algorithm is slow convergence.This can be ameliorated by Nesterov acceleration (Nesterov, 2013).There is also the question of how fast one should send ρ to ∞.Although no optimal schedule is known, simple numerical experiments usually yield a good choice.Finally, soft constraints can be achieved by stopping the steady increase of ρ at a finite value.

Proposed Framework
This simple version of distance majorization can be generalized in various ways.For instance, it can be expanded to multiple constraint sets.In practice, at most two constraint sets usually suffice.Another generalization is to replace the constraint x ∈ S by the constraint Dx ∈ S, where D is a compatible matrix.Again, the original case D = I is allowed.By analogy with the fused lasso of Tibshirani et al. (2005), we will call the matrix D a fusion matrix.This paper is devoted to the study of the general problem of minimizing a differentiable function f (x) subject to r fused constraints D i x ∈ S i .We will approach this problem by extending the proximal distance method.For a fixed penalty constant ρ, the objective function and its MM surrogate now become where P i (y) denotes the projection of y onto S i .Any or all of the fusion matrices D i can be the identity I.
Fortunately, we can simplify the problem by defining S to be the Cartesian product r i=1 S i and D to be the stacked matrix Our objective and surrogate then revert to the less complicated forms (2) where P(x) is the Cartesian product of the projections P i (x).Note that all closed sets S i with simple projections, including sparsity sets, are fair game.

Contributions
In the framework described above, we summarize the contributions of the current paper.
(a) Section 2 describes different solution algorithms for minimizing the penalized loss h ρ (x).
Our first algorithm is based on Newton's method applied to the surrogate g ρ (x | x n ).
For some important problems, Newton's method reduces to least squares.Our second method is a steepest descent algorithm on g ρ (x | x n ) tailored to high dimensional problems.
(b) For a sufficiently large ρ, we show that when f (x) and S are convex and f (x) possesses a unique minimum point y ∈ D −1 (S), the penalized loss h ρ (x) attains its minimum value.This is the content of Proposition 3.1.Similarly, Proposition 3.2 shows that the surrogate g ρ (x | x n ) also attains its minimum.
(c) If in addition f (x) is differentiable, then Proposition 3.3 demonstrates that the MM iterates x n for minimizing h ρ (x) satisfy where z ρ minimizes h ρ (x).If f (x) is also L-smooth and µ-strongly convex, then Proposition 3.4 shows that z ρ is unique and the iterates x n converge to z ρ at a linear rate.
(d) More generally, Proposition 4.2 shows that the iterates x n of a generic MM algorithm for minimizing a coercive subanalytic function h(x) with a good surrogate converge to a stationary point.Our objectives and their surrogates fall into this category.
(e) Finally, we discuss a competing alternating direction method of multipliers (ADMM) algorithm and note its constituent updates.Our extensive numerical experiments compare the two proximal distance algorithms to ADMM.We find that proximal distance algorithms are competitive with and often superior to ADMM in terms of accuracy and running time.

Different Solution Algorithms
Unless f (x) is a convex quadratic, exact minimization of the surrogate g ρ (x | x n ) is likely infeasible.As we have already mentioned, to reduce the objective h ρ (x) in (2), it suffices to reduce the surrogate (3).For the latter task, we recommend Newton's method on small and intermediate-sized problems and steepest descent on large problems.The exact nature of these generic methods are problem dependent.The following section provides a high-level overview of each strategy and we defer details on our later numerical experiments to the appendices.

Newton's Method and Least Squares
Unfortunately, the proximal operator prox ρ −1 f (y) is no longer relevant in calculating the MM update x n+1 .When f (x) is smooth, Newton's method for the surrogate g ρ (x | x n ) employs the update where is the Hessian.To enforce the descent property, it is often prudent to substitute a positive definite approximation H n for d 2 f (x n ).In statistical applications, the expected information matrix is a natural substitute.It is also crucial to retain as much curvature information on f (x) as possible.Newton's method has two drawbacks.First, it is necessary to compute and store d 2 f (x n ).This is mitigated in statistical applications by the substitution just mentioned.Second, there is the necessity of solving a large linear system.Fortunately, the matrix H n + ρD t D is often well-conditioned, for example, when D has full column rank and D t D is positive definite.The method of conjugate gradients can be called on to solve the linear system in this ideal circumstance.
To reduce the condition number of the matrix H n + ρD t D even further, one can sometimes rephrase the Newton step as iteratively reweighted least squares.For instance, in a generalized linear model, the gradient ∇f (x) and the expected information H can be written as where r is a vector of standardized residuals, Z is a design matrix, and W is a diagonal matrix of case weights (Lange, 2010;Nelder and Wedderburn, 1972).The Newton step is now equivalent to minimizing the least squares criterion In this context a version of the conjugate gradient algorithm adapted to least squares is attractive.The algorithms LSQR (Paige and Saunders, 1982) and LSMR (Fong and Saunders, 2011) perform well when the design is sparse or ill conditioning is an issue.

Proximal Distance by Steepest Descent
In high-dimensional optimization problems, gradient descent is typically employed to avoid matrix inversion.Determination of an appropriate step length is now a primary concern.In the presence of fusion constraints Dx ∈ S and a convex quadratic loss f (x) = 1 2 x t Ax+b t x, the gradient of the proximal distance objective at x n amounts to For the steepest descent update x n+1 = x n − t n v n , one can show that the optimal step length is . This update obeys the descent property and avoids matrix inversion.One can also substitute a local convex quadratic approximation around x n for f (x).If the approximation majorizes f (x), then the descent property is preserved.In the failure of majorization, the safeguard of step halving is trivial to implement.
In addition to Nesterov acceleration, gradient descent can be accelerated by the subspace MM technique (Chouzenoux et al., 2010).Let G n be the matrix with k columns determined by the k most current gradients of the objective h ρ (x), including ∇h ρ (x n ).Generalizing our previous assumption, suppose f (x) has a quadratic surrogate with Hessian H n at x n .Overall we get the quadratic surrogate We now seek the best linear perturbation x n + G n β of x n by minimizing q ρ (x n + G n β | x n ) with respect to the coefficient vector β.To achieve this end, we solve the stationary equation The indicated matrix inverse is just k × k.

ADMM
ADMM (alternating direction method of multipliers) is a natural competitor to the proximal distance algorithms just described (Hong et al., 2016).ADMM is designed to minimize functions of the form f (x) + g(Dx) subject to x ∈ C, where C is closed and convex.Splitting variables leads to the revised objective f (x) + g(y) subject to x ∈ C and y = Dx.

ADMM invokes the augmented Lagrangian
with Lagrange multiplier λ and step length µ > 0. At iteration n + 1 of ADMM one calculates successively (5) Update (4) succumbs to Newton's method when f (x) is smooth and C = R p , and update (5) succumbs to a proximal map of g(y).Update (6) of the Lagrange multiplier λ amounts to steepest ascent on the dual function.A standard extension to the scheme in (4) through ( 6) is to vary the step length µ by considering the magnitude of residuals (Boyd et al., 2011).For example, letting r n = Dx − y and s n = µD t (y n−1 − y n ) denote primal and dual residuals at iteration n, we make use of the heuristic which (a) keeps the primal and dual residuals within an order of magnitude of each other, (b) makes ADMM less sensitive to the choice of step length, and (c) improves convergence.
Our problem conforms to the ADMM paradigm when S is equal to the Cartesian product r i=1 S i and g(y) = ρ 2 dist(y, S) 2 .Fortunately, the y update (5) reduces to a simple formula (Bauschke and Combettes, 2017).To derive this formula, note that the proximal map y = prox αg (z) satisfies the stationary condition for any z, including z = Dx n+1 + λ n , and any α, including α = ρ/µ.Since the projection map P(y) has the constant value P(z) on the line segment [z, P(z)], the value satisfies the stationary condition.Because the explicit update (5) for y decreases the Lagrangian even when S is nonconvex, we will employ it generally.The x update (4) is given by the proximal map prox µ −1 f (λ n − y n ) when S = R p and D = I.Otherwise, the update of x is more problematic.Assuming f (x) is smooth and S = R p , Newton's method gives the approximate update Our earlier suggestion of replacing d 2 f (x n ) by a positive definite approximation also applies here.Let us emphasize that ADMM eliminates the need for distance majorization.Although distance majorization is convenient, it is not necessarily a tight majorization.Thus, one can hope to see gains in rates of convergence.Balanced against this positive is the fact that ADMM is often slow to converge to high accuracy.

Proximal Distance Iteration
We conclude this section by describing proximal distance algorithms in pseudocode.As our theoretical results will soon illustrate, the choice of penalty parameter ρ is tied to the convergence rate of any proximal distance algorithm.Unfortunately, a large value for ρ is necessary for the iterates to converge to the constraint set S. We ameliorate this issue by slowly sending ρ → ∞ according to annealing schedules from the family of geometric progressions ρ(t) = r t−1 with t ≥ 1.Here we parameterize the family by an initial value ρ = 1 and a multiplier r > 1.Thus, our methods approximate solutions to min f (x) subject to Dx ∈ S by solving a sequence of increasingly penalized subproblems, min x h ρ(t) (x).In practice we can only solve a finite number of subproblems so we prescribe the following convergence criteria Condition ( 7) is a guarantee that a solution estimate x n is close to a stationary point after n inner iterations for the fixed value of ρ = ρ(t).In conditions ( 8) and (9), the vector x t denotes the δ h -optimal solution estimate once condition (7) is satisfied for a particular subproblem along the annealing path.Condition (8) is a guarantee that solutions along the annealing path adhere to the fusion constraints at level δ d .In general, condition (8) can only be satisfied for large values of ρ(t).Finally, condition (9) is used to terminate the annealing process if the relative progress made in decreasing the distance penalty becomes too small as measured by δ q .Algorithm 1 summarizes the flow of proximal distance iteration, which uses Nesterov acceleration in inner iterations.Warm starts in solving subsequent subproblem are implicit in our formulation.

Convergence Analysis: Convex Case
Let us begin by establishing the existence of a minimum point.Further constraints on x beyond those imposed in the distance penalties are ignored or rolled into the essential domain of f (x) when f (x) is convex.As noted earlier, we can assume a single fusion matrix D and a single closed convex constraint set S. In such setting we have the following result.Proofs are deferred to Section 7. Proposition 3.1 Suppose the convex function f (x) on R p possesses a unique minimum point y on the closed convex set T = D −1 (S).Then for all sufficiently large ρ, the objective h ρ (x) = f (x) + ρ 2 dist(Dx, S) 2 is coercive and therefore attains its minimum value.Next we show that the majorization surrogate defined in (3) attains its minimum value for large enough ρ.Proposition 3.2 Under the conditions of Proposition 3.1, for sufficiently large ρ, every surrogate g ρ (x | x n ) = f (x) + ρ 2 Dx − P(Dx n ) 2 is coercive and therefore attains its minimum value.If for all x and some positive semidefinite matrix A and subgradient v n at x n , and if the inequality u t Au > 0 holds whenever Du = 0 and u = 0, then for ρ sufficiently large, g ρ (x | x n ) is strongly convex and hence coercive.
We continue to assume that f (x) and S are convex and that a minimum point of the surrogate g ρ (x | x n ) is available.Uniqueness of x n+1 holds when g ρ (x | x n ) is strictly convex.The constraint set S is implicitly captured by the essential domain of f (x).
Our earlier research shows that moving some constraints to the essential domain of f (x) is sometimes helpful (Keys et al., 2019;Lange, 2016).In any event, in our ideal convex setting we have a first convergence result for fixed ρ.
Proposition 3.3 Supposes a) that S is closed and convex, b) that the loss f (x) is convex and differentiable, and c) that the constrained problem possesses a unique minimum point.For ρ sufficiently large, let z ρ denote a minimal point of the objective h ρ (x) defined by equation ( 2).Then the MM iterates (10) satisfy Furthermore, the iterate values h ρ (x n ) systematically decrease.
In even more restricted circumstances, one can prove linear convergence of function values in the framework of (Karimi et al., 2016).
Proposition 3.4 Suppose that S is a closed and convex set and that the loss f (x) is Lsmooth and µ-strongly convex.Then the objective h ρ (x) = f (x) + ρ 2 dist(Dx, S) 2 possesses a unique minimum point z ρ , and the proximal distance iterates x n satisfy Convergence of ADMM is well studied in the optimization literature (Beck, 2017).Appendix A summarizes the main findings.

Convergence Analysis: General Case
We now depart the comfortable confines of convexity.Let us first review the notion of a Fréchet subdifferential (Kruger, 2003) For a locally Lipschitz and directionally differentiable function, the Fréchet subdifferential becomes Here d u h(x) is the directional derivative of h(x) at x in the direction u.This result makes it clear that at a critical point, all directional derivatives are flat or point uphill.
We will also need some notions from algebraic geometry (Bochnak et al., 2013).For simplicity we focus on the class of semialgebraic functions and the corresponding class of semialgebraic subsets of R p .The latter is the smallest class that: (a) contains all sets of the form {x : q(x) > 0} for a polynomial q(x) in p variables, (b) is closed under the formation of finite unions, finite intersections, set complementation, and Cartesian products.
A function a : R p → R r is said to be semialgebraic if its graph is a semialgebraic set of R p+r .The class of real-valued semialgebraic functions contains all polynomials p(x) and all 0/1 indicators of algebraic sets.It is closed under the formation of sums and products and therefore constitutes a commutative ring with identity.The class is also closed under the formation of absolute values, reciprocals when a(x) = 0, nth roots when a(x) ≥ 0, and maxima max{a(x), b(x)} and minima min{a(x), b(x)}.Finally, the composition of two semialgebraic functions is semialgebraic.
For our purposes it is crucial that the Euclidean distance dist(x, S) to a semialgebraic set S is a semialgebraic function.In view of the closure properties of such functions, the function ρ 2 dist(Dx, S) 2 is also semialgebraic.Sets such as the nonnegative orthant R p + and the unit sphere S p−1 are semialgebraic.The next proposition buttresses several of our numerical examples.
Proposition 4.1 The order statistics of a finite set {f i (x)} n i=1 of semialgebraic functions are semialgebraic.Hence, sparsity sets are semialgebraic.
The next proposition is an elaboration and expansion of known results (Attouch et al., 2010;Bolte et al., 2007;Cui et al., 2018;Kang et al., 2015;Le Thi et al., 2018) and was featured in our previous paper (Keys et al., 2019).Proposition 4.2 In an MM algorithm suppose the objective h(x) is coercive, continuous, and subanalytic and all surrogates g(x | x n ) are continuous, µ-strongly convex, and satisfy the Lipschitz condition Proposition 4.2 applies to proximal distance algorithms under the right hypotheses.Note that semialgebraic sets and functions are automatically subanalytic.Before stating a precise result, let us clarify the nature of the Fréchet subdifferential in the current setting.This entity is determined by the identity Danskin's theorem yields the directional derivative where S(x) is the solution set where the minimum is attained.The Frechet differential holds owing to Corollary 1.12.2 and Proposition 1.17 of (Kruger, 2003) since dist(Dx, S) 2 is locally Lipshitz.The latter fact follows from the identity ) is Lipschitz and bounded on bounded sets.
In any event, a stationary point x satisfies 0 = ∇f (x) + ρD t (Dx − z) for all z ∈ S(x).As we expect, the stationary condition is necessary for x to furnish a global minimum.Indeed, if it fails, take z ∈ S(x) with surrogate satisfying ∇g ρ (x | x) = 0. Then the negative gradient −∇g ρ (x | x) is a descent direction for g ρ (x | x), which majorizes h ρ (x).Hence, −∇g ρ (x | x) is also a descent direction for h ρ (x).This conclusion is inconsistent with x being a local minimum of the objective.
The next proposition proves convergence for a wide class of fused models.
Proposition 4.3 Suppose in our proximal distance setting that ρ is sufficiently large, the closed constraint sets S i and the loss f (x) are semialgebraic, and f (x) is differentiable with a locally Lipschitz gradient.Under the coercive assumption made in Proposition 3.2, the proximal distance iterates x n converge to a stationary point of the objective h ρ (x).
For sparsity constrained problems, one can establish a linear rate of convergence under the right hypotheses.
Proposition 4.4 Suppose in our proximal distance setting that ρ is sufficiently large, the constraint set is a sparsity set S k , and the loss f (x) is semialgebraic, strongly convex, and possesses a Lipschitz gradient.Then the proximal distance iterates x n converge linearly to a stationary point x ∞ provided Dx ∞ has k unambiguous largest components in magnitude.When the rows of D are unique, the complementary set of points x where Dx has k ambiguous largest components in magnitude has Lebesgue measure 0.

Numerical Examples
This section considers five concrete examples of constrained optimization amenable to distance majorization with fusion constraints, with D denoting the fusion matrix in each problem.In each case, the loss function is both strongly convex and differentiable.The specific examples that we consider are the metric projection problem, convex regression, convex clustering, image denoising with a total variation penalty, and projection of a matrix to one with a better condition number.Each example is notable for the large number of fusion constraints and projections to convex constraint sets, except in convex clustering.In convex clustering we encounter a sparsity constraint set.Quadratic loss models feature prominently in our examples.Interested readers may consult our previous work for nonconvex examples with D = I (Keys et al., 2019;Xu et al., 2017).

Mathematical Descriptions
Here we provide the mathematical details for each example.

Metric Projection
Solutions to the metric projection problem restore transitivity to noisy distance data for the m nodes of a graph (Brickell et al., 2008;Sra et al., 2005).The data are encoded in an m×m dissimilarity matrix Y = (y ij ) with nonnegative weights in the matrix W = (w ij ).The metric projection problem requires finding the symmetric semi-metric X = (x ij ) minimizing subject to all m 2 nonnegativity constraints x ij ≥ 0 and all 3 m 3 triangle inequality constraints x ij − x ik − x kj ≤ 0. The diagonal entries of Y , W , and X are zero by definition.The fusion matrix D has m 2 + 3 m 3 rows, and the projected value of DX must fall in the set S of symmetric matrices satisfying all pertinent constraints.
One can simplify the required projection by stacking the nonredundant entries along each successive column of X to create a vector x with m 2 entries.This captures the lower triangle of X.The sparse matrix D is correspondingly redefined to be [ m 2 + 3 m 3 ] × m 2 .These maneuvers simplfy constraints to Dx ≥ 0, and projection involves sending each entry u of Dx to max{0, u}.Putting everything together, the objective to minimize is where D consists of blocks T and I p 2 and p 1 and p 2 count the number of triangle inequality and nonnegativity constraints, respectively.The linear system (I + ρD t D)x = b appears in both the MM and ADMM updates for x n .Application of the Woodbury and Sherman-Morrison formulas yield an exact solution to the linear system and allow one to forgo iterative methods.The interested reader may consult Appendix B for further details.

Convex Regression
Convex regression is a nonparametric method for estimating a regression function under shape constraints.Given m responses y i and corresponding predictors x i ∈ R d , the goal is to find the convex function ψ(x) minimizing the sum of squares . Asymptotic and finite sample properties of this convex estimator have been described in detail by Seijo and Sen (2011).The convex regression program can be restated as the finite dimensional problem of finding the value θ i and subgradient ξ i ∈ R d of ψ(x) at each sample point (y i , x i ).Convexity imposes the supporting hyperplane constraint θ j +ξ t j (x i −x j ) ≤ θ i for each pair i = j.Thus, the problem becomes one of minimizing 1 2 y−θ 2 subject to these m(m − 1) inequality constraints.In the proximal distance framework, we must minimize where D = [A B] encodes the required fusion matrix.The reader may consult Appendix C for a description of each algorithm map.

Convex Clustering
Convex clustering of m samples based on d features can be formulated in terms of the regularized objective based on columns x i and u i of X ∈ R d×m and U ∈ R d×m , respectively.Here each x i ∈ R d is a sample feature vector and the corresponding u i represents its centroid assignment.
The predetermined weights w ij have a graphical interpretation under which similar samples have positive edge weights w ij and distant samples have 0 edge weights.The edge weights are chosen by the user to guide the clustering process.In general, minimization of F γ (U ) separates over the connected components of the graph.To allow all sample points to coalesce into a single cluster, we assume that the underlying graph is connected.The regularization parameter γ > 0 tunes the number of clusters in a nonlinear fashion and potentially captures hierarchical information.Previous work establishes that the solution path U (γ) varies continuously with respect to γ (Chi and Lange, 2015).Unfortunately, there is no explicit way to determine the number of clusters entailed by a particular value of γ.
Alternatively, we can attack the problem using sparsity and distance majorization.Consider the penalized objective The fusion matrix D has m 2 columns w ij (e i − e j ) and serves to map the centroid matrix U to a d × m 2 matrix V encoding the weighted differences w ij (u i − u j ).The members of the sparsity set S k are d × m 2 matrices with at most k non-zero columns.Projection of U D onto the closed set S k forces some centroid assignments to coalesce, and is straightforward to implement by sorting the Euclidean lengths of the columns of U D and sending to 0 all but the k most dominant columns.Ties are broken arbitrarily.
Our sparsity-based method trades the continuous penalty parameter γ > 0 in the previous formulation for an integer sparsity index k ∈ {0, 1, 2, . . ., m 2 }.For example with k = 0, all differences u i − u j are coerced to 0, and all sample points cluster together.The other extreme k = m 2 assigns each point to its own cluster.The size of the matrices D and U D can be reduced by discarding column pairs corresponding to 0 weights.Appendix D describes the projection onto sparsity sets and provides further details.

Total Variation Image Denoising
To approximate an image U from a noisy input W matrix, Rudin et al. (1992) regularize a loss function f (U ) by a total variation (TV) penalty.After discretizing the problem, the least squares loss leads to the objective where U , W ∈ R m×p are rectangular monochromatic images and γ controls the strength of regularization.The anisotropic norm is often preferred because it induces sparsity in the differences.Here D p is the forward difference operator on p data points.Stacking the columns of U into a vector u = vec(U ) allows one to identify a fusion matrix D and write TV 1 (U ) compactly as TV 1 (u) = Du 1 .In this context we reformulate the denoising problem as minimizing f (u) subject to the set constraint Du 1 ≤ γ.This revised formulation directly quantifies the quality of a solution in terms of its total variation and brings into play fast pivot-based algorithms for projecting onto multiples of the 1 unit ball (Condat, 2016).Appendix E provides descriptions of each algorithm.

Projection of a Matrix to a Good Condition Number
Consider an m×p matrix M with m ≥ p and full singular value decomposition M = U ΣV t .The condition number of M is the ratio σ max /σ min of the largest to the smallest singular value of M .We denote the diagonal of Σ as σ.Owing to the von Neumann-Fan inequality, the closest matrix N to M in the Frobenius norm has the singular value decomposition N = U XV t , where the diagonal x of X satisfies inequalities pertinent to a decent condition number (Borwein and Lewis, 2010).Suppose c ≥ 1 is the maximum condition number.Then every pair (x i , x j ) satisfies x i − cx j ≤ 0. Note that x i − cx i > 0 if and only if x i < 0. Thus, nonnegativity of the entries of x is enforced.The proximal distance approach to the condition number projection problem invokes the objective and majorization at iteration n, where q nij = min{x ni − cx nj , 0}.We can write the majorization more concisely as where vec Q n stacks the columns of Q n = (q nij ) and the p 2 × p fusion matrix D satisfies (Dx) k = x i − cx j for each component k.The minimum of the surrogate occurs at the point x n+1 = (A t ρ A ρ ) −1 A t ρ r n .This linear system can be solved exactly.Appendix F provides additional details.

Numerical Results
Our numerical experiments compare various strategies for implementing Algorithm 1.We consider two variants of proximal distance algorithms.The first directly minimizes the majorizing surrogate (MM), while the second performs steepest descent (SD) on it.In addition to the aforementioned methods, we tried the subspace MM algorithm described in Section 2.2.Unfortunately, this method was outperformed in both time and accuracy comparisons by Nesterov accelerated MM; the MM subspace results are therefore omitted.We also compare our proximal distance approach to ADMM as described in Section 2.3.In many cases updates require solving a large linear system; we found that the method of conjugate gradients sacrificed little accuracy and largely outperformed LSQR and therefore omit comparisons.The clustering and denoising examples are exceptional in that the associated matrices D t D are sufficiently ill-conditioned to cause failures in conjugate gradients.Table 1 summarizes choices in control parameters across each example.
We now explain example by example the implementation details behind our efforts to benchmark the three strategies (MM, SD, and ADMM) in implementing Algorithm 1.In each case we initialize the algorithm with the solution of the corresponding unconstrained problem.Performance is assessed in terms of speed in seconds or milliseconds, number of iterations until convergence, the converged value of the loss f (x), and the converged distance to the constraint set dist(Dx, S), as described in Algorithm 1.Additional metrics are highlighted where applicable.The term inner iterations refers to the number of iterations to solve a penalized subproblem argmin h ρ (x) for a given ρ whereas outer iterations count the total number of subproblems solved.Lastly, we remind readers that the approximate solution to argmin h ρ(t) (x) is used as a warm start in solving argmin h ρ(t+1) (x).

Metric Projection.
In our comparisons, we use input matrices Y ∈ R m×m whose iid entries y ij are drawn uniformly from the interval [0, 10] and set weights w ij = 1.Each algorithm is allotted a maximum of 200 outer and 10 5 inner iterations, respectively, to achieve a gradient norm of δ h = 10 −3 and distance to feasibility of δ d = 10 −2 .The relative change parameter is set to δ q = 0 and the annealing schedule is set to ρ(t) = min{10 8 , 1.2 t−1 } for the proximal distance methods.Table 2 summarizes the performance of the three algorithms.Best values appear in boldface.All three algorithms converge to a similar solution as indicated by the final loss y − x 2 and distance values.It is clear that SD matches or outperforms MM and ADMM approaches on this example.Notably, the linear system appearing in the MM update admits an exact solution (see Appendix B.5), yet SD has a faster time to  solution with fewer iterations taken.The selected convergence metrics in Figure 1 vividly illustrate stability of solutions x ρ = argmin h ρ (x) along an annealing path from ρ = 1 to ρ = 1.2 40 ≈ 1470.Specifically, solving each penalized subproblem along the sequence results in marginal increase in the loss term with appreciable decrease in the distance penalty.
Except for the first outer iteration, there is minimal decrease of the loss, distance penalty, or penalized objective within a given outer iteration even as the gradient norm vanishes.
The observed tradeoff between minimizing a loss model and minimizing a nonnegative penalty is well-known in penalized optimization literature (Beltrami, 1970; Lange, 2016, see Proposition 7.6.1 on p. 183).
Figure 1: Loss, distance, penalized objective, and gradient norm for SD on metric projection problem with 32 nodes, labeled by outer iteration.

Convex Regression.
In our numerical examples the observed functional values y i are independent Gaussian deviates with means ψ(x i ) and common variance σ 2 = 0.1.The predictors are iid deviates sampled from the uniform distribution on [−1, 1] d .We choose the simple convex function ψ(x i ) = x i 2 for our benchmarks for ease in interpretation; the interested reader may consult the work of Mazumder et al. (2019) for a detailed account of the applicability of the technique in general.Each algorithm is allotted a maximum of 200 outer and 10 4 inner iterations, respectively, to converge with δ h = 10 −3 , δ d = 10 −2 , and δ q = 10 −6 .The annealing schedule is set to ρ(t) = min{10 8 , 1.2 t−1 }.
Table 3 demonstrates that although the SD approach is appreciably faster than both MM and ADMM, the latter appear to converge on solutions with marginal improvements in minimizing the loss y − θ 2 , distance, and mean squared error (MSE) measured using ground truth functional values ψ(x i ) and estimates θ i .Interestingly, increasing both the number of features and samples does not necessarily increase the amount of required computational time in using a proximal distance approach; for example, see results with d = 2 and d = 20 features.This may be explained by sensitivity to the annealing schedule.

Convex Clustering.
To evaluate the performance of the different methods on convex clustering, we consider a mixture of simulated data and discriminant analysis data from the UCI Machine Learning Repository (Dua and Graff, 2019).The simulated data in gaussian300 consists of 3 Gaussian clusters generated from bivariate normal distributions with means µ = (0.0, 0.0) t , (2.0, 2.0) t , and (1.8, 0.5) t , standard deviation σ = 0.1, and class sizes n 1 = 150, n 2 = 50, n 3 = 100.This easy dataset is included to validate Algorithm 2 described later as a reasonable solution path heuristic.The data in iris and zoo are representative of clustering with purely continuous or purely discrete data, respectively.In these two datasets, samples with same class label form a cluster.Finally, the simulated data spiral500 is a classic example that thwarts k-means clustering.Each algorithm is allotted a maximum of 10 4 inner iterations to solve a ρ-penalized subproblem at level δ h = 10 −2 .The annealing schedule is set to ρ(t) = min{10 8 , 1.2 t−1 } over 100 outer iterations with δ d = 10 −5 and δ q = 10 −6 .
Because the number of clusters is usually unknown, we implement the search heuristic outlined in Algorithm 2. The idea behind the heuristic is to gradually coerce clustering without exploring the full range of the hyperparameter k.As one decreases the number of admissible nonzero centroid differences k from k max to 0, sparsity (1−k/k max ) in the columns of U D increases to reflect coalescing centroid assignments.Thus, Algorithm 1 generates a list of candidate clusters that can be evaluated by various measures of similarity (Vinh et al., 2010).For example, the adjusted Rand index (ARI) provides a reasonable measure of the distance to the ground truth in our examples as it accounts for both the number of identified clusters and cluster assignments.We also report the related normed Mutual Information (NMI).The ARI takes values on [−1, 1] whereas NMI appears on a [0, 1] scale.
ADMM, as implemented here, is not remotely competitive on these examples given its extremely long compute times and failure to converge in some instances.These times are only exacerbated by the search heuristic and therefore omit ADMM from this example.The findings reported in Table 4 indicate the same accuracy for MM (using LSQR) and SD as measured by loss and distance to feasibility.Here we see that the combination of the  proximal distance algorithms and the overall search heuristic (Algorithm 2) yields perfect clusters in the gaussian300 example on the basis of ARI and NMI.To its disadvantage, the search heuristic is greedy and generally requires tuning.Both MM and SD achieve similar clusterings as indicated by ARI and NMI.Notably, SD generates candidate clusterings faster than MM.

Total Variation Image Denoising.
To evaluate our denoising algorithm, we consider two standard test images, cameraman and peppers gray.White noise with σ = 0.2 is applied to an image and then reconstructed using our proximal distance algorithms.Only MM and SD are tested with a maximum of 100 outer and 10 4 inner iterations, respectively, and convergence thresholds δ h = 10 −1 , δ d = 10 −1 , and δ q = 10 −6 .A moderate schedule ρ(t) = {10 8 , 1.5 t−1 } performs well even with such lax convergence criteria.Table 5 reports convergence metrics and image quality indices, MSE and Peak Signal-to-Noise Ratio (PSNR).Timings reflect the total time spent generating solutions, starting from a 0% reduction in the total variation of the input image U up to 90% reduction in increments of 10%.Explicitly, we take γ 0 = TV 1 (U ) and vary the control parameter γ = (1 − s)γ 0 with s ∈ [0, 1] to control the strength of denoising.Figure 5.2.4 depicts the original and reconstructed images along the solution path.Here c(M ) is the condition number of the input matrix, a is the decrease factor, and c(X) is the condition number of the solution.

Projection of a Matrix to a Good Condition Number.
We generate base matrices M ∈ R p×p as random correlation matrices using MatrixDepot.jl(Zhang and Higham, 2016), which relies on Davies' and Higham's refinement (Davies and Higham, 2000) of the Bendel-Mickey algorithm (Bendel and Mickey, 1978).Simulations generate matrices with condition numbers c(M ) in the set {119, 1920, 690}.Our subsequent analyses target condition number decreases by a factor a such a ∈ {2, 4, 16, 32}.Each algorithm is allotted a maximum of 200 outer and 10 4 inner iterations, respectively with choices δ h = 10 −3 , δ d = 10 −2 , δ q = 10 −6 , and ρ(t) = 1.2 t−1 .Table 6 summarizes the performance of the three algorithms.The quality of approximate solutions is similar across MM, SD, and ADMM in terms of loss, distance, and final condition number metrics.Interestingly, the MM approach requires less time to deliver solutions of comparable quality to SD solutions as the size p of the input matrix M ∈ R p×p increases.

Discussion
We now recapitulate the main findings of our numerical experiments.Tables 2 through  6 show a consistent pattern of superior speed by the steepest descent (SD) version of the proximal distance algorithm.This is hardly surprising since unlike ADMM and MM, SD avoids solving a linear system at each iteration.SD's speed advantage tends to persist even when the linear system can be solved exactly.The condition number example summarized in Table 6 is an exception to this rule.Here the MM updates leverage a very simple matrix inverse.MM is usually faster than ADMM.We attribute MM's performance edge to the extra matrix-vector multiplications involving the fusion matrix D required by ADMM.In fairness, ADMM closes the speed gap and matches MM on convex regression.The choice of annealing schedule can strongly impact the quality of solutions.Intuitively, driving the gradient norm ∇h ρ (x) to nearly 0 for a given ρ keeps the proximal distance methods on the correct annealing path and yields better approximations.Provided the penalized objective is sufficiently smooth, one expects the solution x ρ ≈ argmin h ρ (x) to be close to the solution x ρ ≈ argmin h ρ (x) when the ratio ρ /ρ > 1 is not too large.Thus, choosing a conservative δ h for the convergence criterion ∇h ρ (x) ≤ δ h may guard against a poorly specified annealing schedule.Quantifying sensitivity of intermediate solutions x ρ with respect to ρ is key in avoiding an increase in inner iterations per subproblem; for example, as observed in Figure 1.Given the success of our practical annealing recommendation to overcome the unfortunate coefficients in Propositions 3.3 and 3.4, this topic merits further consideration in future work.
In practice, it is sometimes unnecessary to impose strict convergence criteria on the proximal distance iterates.It is apparent that the convergence criteria on convex clustering and image denoising are quite lax compared to choices in other examples, specifically in terms of δ h .Figure 1 suggests that most of the work in metric projection involves driving the distance penalty downhill rather than in fitting the loss.Surpisingly, Table 4 shows that our strict distance criterion δ d = 10 −5 in clustering is achieved.This implies dist(Dx, S) 2 ≤ 10 −10 on the selected solutions with ρ ≤ 10 8 , yet we only required ∇h ρ (x) ≤ 10 −2 on each subproblem.Indeed, not every problem may benefit from precise solution estimates.The image processing example underscores this point as we are able to recover denoised images with the choices δ h = δ d = 10 −1 .Problems where patterns or structure in solutions are of primary interest may stand to benefit from relaxed convergence criteria.
Our proximal distance method, as described in Algorithm 1, enjoys several advantages.First, fusion constraints fit naturally in the proximal distance framework.Second, proximal distances enjoy the descent property.Third, there is a nearly optimal step size for gradient descent when second-order information is available on the loss.Fourth, proximal distance algorithms are competitive if not superior to ADMM on many problems.Fifth, proximal distance algorithms like iterative hard thresholding rely on set projection and are therefore helpful in dealing with hard sparsity constraints.The main disadvantages of the proximal distance methods are (a) the overall slow convergence due to the loss of curvature information on the distance penalty and (b) the need for a reasonable annealing schedule.In practice, a little experimentation can yield a reasonable schedule for an entire class of problems.Many competing methods are only capable of dealing with soft constraints imposed by the lasso and other convex penalties.To their detriment, soft constraints often entail severe parameter shrinkage and lead to an excess of false positives in model selection.
Throughout this manuscript we have stressed the flexibility of the proximal distance framework in dealing with a wide range of constraints as a major strength.From our point of view, proximal distance iteration adequately approximates feasible, locally optimal solutions to constrained optimization problems for well behaved constraint sets, for instance convex sets or semialgebraic sets.Combinatorially complex constraints or erratic loss functions can cause difficulties for the proximal distance methods.The quadratic distance penalty dist(Dx, S) 2 is usually not an issue, and projection onto the constraint should be fast.Poor loss functions may either lack second derivatives or may possess a prohibitively expensive and potentially ill-conditioned Hessian d 2 f (x).In this setting techniques such as coordinate descent and regularized and quasi-Newton methods are viable alternatives for minimizing the surrogate g ρ (x | x n ) generated by distance majorization.In any event, it is crucial to design a surrogate that renders each subproblem along the annealing path easy to solve.This may entail applying additional majorizations in f (x).Balanced against this possibility is the sacrifice of curvature information with each additional majorization.
We readily acknowledge that other algorithms may perform better than MM and proximal distance algorithms on specific problems.The triangle fixing algorithm for metric projection is a case in point (Brickell et al., 2008), as are the numerous denoising algorithms based on the 1 norm.This objection obscures the generic utility of the proximal distance principle.ADMM can certainly be beat on many specific problems, but nobody seriously suggests that it be rejected across the board.Optimization, particularly constrained optimization, is a fragmented subject, with no clear winner across problem domains.Generic methods serve as workhorses, benchmarks, and backstops.
As an aside, let us briefly note that ADMM can be motivated by the MM principle, which is the same idea driving proximal distance algorithms.The optimal pair (x, y) and λ furnishes a stationary point of the Lagrangian.Because the Lagrangian is linear in λ, its maximum for fixed (x, y) is ∞.To correct this defect, one can add a viscosity minorization to the Lagrangian.This produces the modified Lagrangian The penalty term has no impact on the x and y updates.However, the MM update for λ is determined by the stationary condition so that The choice α = 1 gives the standard ADMM update.Thus, the ADMM algorithm alternates decreasing and increasing the Lagrangian in a search for the saddlepoint represented by the optimal trio (x, y, λ).
In closing we would like to draw the reader's attention to some generalizations of the MM principle and connections to other well-studied algorithm classes.For instance, a linear fusion constraint Dx ∈ S can in principle by replaced by a nonlinear fusion constraint M (x) ∈ S. The objective and majorizer are then The objective has gradient g = ∇f (x)+ρdM (x) t {M (x)−P S [M (x)]}.The second differential of the majorizer is approximately d 2 f (x)+ρdM (x) t dM (x) for M (x) close to P S [M (x)].Thus, gradient descent can be implemented with step size assuming the denominator is positive.Algebraic penalties such as g(x) 2 reduce to distance penalties with constraint set {0}.The corresponding projection operator sends any vector y to 0 so that the algebraic penalty g(x) 2 = dist[g(x), {0}] 2 .This observation is pertinent to constrained least squares with g(x) = d − Cx (Golub and Van Loan, 1996).The proximal distance surrogate can be expressed as and minimized by standard least squares algorithms.No annealing is necessary.Inequality constraints g(x) ≤ 0 behave somewhat differently.The proximal distance majorization dist 2 is not the same as the Beltrami quadratic penalty g(x) 2 + (Beltrami, 1970).However, the standard majorization (Lange, 2016)

Proofs
In this section we provide proofs for the convergence results discussed in Section 3 and Section 4 for the convex and noncovex cases, respectively.

Proposition 3.1
Proof Without loss of generality we can translate the coordinates so that y = 0. Let B be the unit sphere {x : x = 1}.Our first aim is to show that h ρ (x) > f (0) throughout B. Consider the set B ∩ T , which is possibly empty.On this set the infimum b of f (x) is attained, so b > f (0) by assumption.The set B\T will be divided into two regions, a narrow zone adjacent to T and the remainder.Now let us show that there exists a δ > 0 such that h ρ (x) ≥ f (x) ≥ f (0) + δ for all x ∈ B with dist(Dx, S) ≤ δ.If this is not so, then there exists a sequence x n ∈ B with f (x n ) < f (0) + 1 k and dist(Dx n , S) ≤ 1 k .By compactness, some subsequence of x n converges to z ∈ B ∩ T with f (z) ≤ f (0), contradicting the uniqueness of y.Finally, let a = min x∈B f (x).To deal with the remaining region take ρ large enough so that a + ρ 2 δ 2 > f (0).For such ρ, h ρ (x) > f (0) everywhere on B. It follows that on the unit ball {x : x ≤ 1}, h ρ (x) is minimized at an interior point.Because h ρ (x) is convex, a local minimum is necessarily a global minimum.

Proposition 3.2
Proof The first assertion follows from the bound g ρ (x | x n ) ≥ h ρ (x).To prove the second assertion, we note that it suffices to prove the existence of some constant ρ > 0 such that the matrix A + ρD t D is positive definite (Debreu, 1952).If no choice of ρ renders A + ρD t D positive definite, then there is a sequence of unit vectors u m and a sequence of scalars ρ m tending to ∞ such that By passing to a subsequence if needed, we may assume that the sequence u m converges to a unit vector u.On the one hand, because D t D is positive semidefinite, inequality (11) compels the conclusions u t m Au m ≤ 0, which must carry over to the limit.On the other hand, dividing inequality (11) by ρ m and taking limits imply u t D t Du ≤ 0 and therefore Du = 0.Because the limit vector u violates the condition u t Au > 0, the required ρ > 0 exists.

Proposition 3.3
Proof Systematic decrease of the iterate values h ρ (x n ) is a consequence of the MM principle.The existence of z ρ follows from Proposition 3.1.To prove the stated bound, first observe that the function g ρ (x | x n ) − ρ 2 Dx 2 is convex, being the sum of the convex function f (x) and a linear function.Because ∇g ρ (x n+1 | x n ) t (x − x n+1 ) ≥ 0 for any x in S, the supporting hyperplane inequality implies that or equivalently Now note that the difference has gradient ∇d(x | y) = P(x) − P(y).
Because P(x) is non-expansive, the gradient ∇d(x | y) is Lipschitz with constant 1.The tangency conditions d(y | y) = 0 and ∇d(y | y) = 0 therefore yield for all x.At a minimum z ρ of h ρ (x), combining inequalities ( 12) and ( 13) gives Adding the result over n and invoking the descent property h ρ (x n+1 ) ≤ h ρ (x n ), telescoping produces the desired error bound This is precisely the asserted bound.

Proposition 3.4
Proof The existence and uniqueness of z ρ are obvious.The remainder of the proof hinges on the facts that h ρ (x) is µ-strongly convex and the surrogate g ρ (x | w) is L+ρ D 2 -smooth for all w.The latter assertion follows from These facts together with ∇g ρ (z The strong convexity condition It follows that ∇h ρ (x) ≥ µ 2 x − z ρ .This last inequality and inequality ( 14) produce the Polyak-Lojasiewicz bound the Polyak-Lojasiewicz bound gives Rearranging this inequality yields which can be iterated to give the stated bound.

Proposition 4.1
Proof The first claim is true owing to the inclusion-exclusion formula and the previously stated closure properties.For n = 3 and k = 2 the inclusion-exclusion formula reads f To prove the second claim, note that a sparsity set in R p with at most k nontrivial coordinates can be expressed as the zero set {x : y (p−k) = 0}, where y i = |x i |.Thus, it is semialgebraic.

Proposition 4.3
Proof To validate the subanalytic premise of Proposition 4.2, first note that semialgebraic functions and sets are automatically subanalytic.The penalized loss Note that Proposition 3.2 furnishes a sufficient condition under which the functions f (x) + ρ 2 Dx 2 + a t x achieve their minima.Linear convergence holds under stronger assumptions.
Proposition A.2 Suppose that S is closed and convex, that the loss f (x) is L-smooth and µ-strongly convex, and that the map determined by D is onto.Then the ADMM iterates converge at a linear rate.Giselsson and Boyd (2016) proved Proposition A.2 by operator methods.A range of convergence rates is specified there.

Appendix B. Additional Details for Metric Projection Example
Given a n × n dissimilarity matrix C = (c ij ) with non-negative weights w ij , our goal is to find a semi-metric X = (x ij ).We start by denoting trivec an operation that maps a symmetric matrix X to a vector x, x = trivec(X) (Figure 3).Then we write the metric projection objective as where c = trivec(C).Here T encodes triangle inequalities and the m i count the number of * x 12 x 13 x 14 x 21 * x 23 x 24 x 31 x 32 * x 34 x 41 x 42 x 43 * x 21 x 31 x 41 x 32 x 34 x contraints of each type.The usual distance majorization furnishes a surrogate The notation P(•, S) denotes projection onto a set S. The fusion matrix D = [T ; I] stacks the two operators; the joint projection operates in a block-wise fashion.

B.1 MM
We rewrite the surrogate explicitly as a least squares problem minimizing Ax − b n 2 2 : , where c ≡ y from the main text.Updating the RHS b n in the linear system reduces to evaluating the projection and copy operations.It is worth noting that triangle fixing algorithms that solve the metric nearness problem operate in the same fashion, except they work one triangle a time.That is, each iteration solves n 3 least squares problems compared to 1 in this formulation.A conjugate gradient type of algorithm solves the normal equations directly using A t A, whereas LSQR type methods use only A and A t .

B.2 Steepest Descent
The updates x n+1 = x n − t n ∇h ρ (x n ) admit an exact solution for the line search parameter t n .Recall the generic formula from the main text:

B.3 ADMM
Taking y as the dual variable and λ as scaled multipliers, the updates for each ADMM block are Finally, the Multipliers follow the standard update.

B.4 Properties of the Triangle Inequality Matrix
These results have been documented before and are useful in designing fast subroutines for Dx and D t Dx.Recall that m counts the number of nodes in the problem and p = m 2 is the number of parameters.In this notation D = T I p and D t D = T t T + I p .
Proposition 1 The matrix T has 3 m 3 rows and m 2 columns.Proof Interpret X as the adjacency matrix for a complete directed graph on m nodes without self-edges.When X is symmetric the number of free parameters is therefore m 2 .An oriented 3-cycle is formed by fixing 3 nodes so there are m 3 such cycles.Now fix the orientation of the 3-cycles and note that each triangle encodes 3 metric constraints.The number of constraints is therefore 3 m 3 .
Proof In view of the previous result, the entries T ij encode whether edge j participates in constraint i.We proceed by induction on the number of nodes m.The base case m = 3 involves one triangle and is trivial.Note that a triangle encodes 3 inequalities.Now consider a complete graph on m nodes and suppose the claim holds.Without loss of generality, consider the collection of 3-cycles oriented clockwise and fix an edge j.Adding a node to the graph yields 2m new edges, two for each of the existing m nodes.This action also creates one new triangle for each existing edge.Thus, edge j appears in 3(m − 2) + 3 = 3[(m + 1) − 2] triangle inequality constraints based on the induction hypothesis.
Proof Interpret the inequality x ij ≤ x ik + x kj with i > k > j as the ordered triple x ij , x ik , x kj .The statement is equivalent to counting a(N ) = number of times x ij appears in position 1, and, b(N ) = number of times x ij appears in position 2 or 3, where N denotes the number of constraints.In view of the previous proposition, it is enough to prove a(N ) = m − 2. Note that a(3) = 1, meaning that x ij appears in position 1 exactly once within a given triangle.Given that an edge (i, j) appears in 3(m − 2) constraints, divide this quantity by the number of constraints per triangle to arrive at the stated result.
Proposition 4 The matrix T has full column rank.
Proof It is enough to show that A = T t T is full rank.The first two propositions imply To compute the off-diagonal entries, fix a triangle and note that two edges i and j appear in all three of its constraints of the form x i ≤ x j + x k .There are three possibilities for a given constraint c: −1, if i LHS, j RHS or vice-versa 1, if i and j both appear on RHS 0, if one of i or j is missing.
It follows that a ij = T i , T j = −1, if edges i and j overlap in constraints 0, otherwise.
By Proposition B.2, an edge i appears in 3(m − 2) constraints.Imposing the condition that edge j also appears reduces this number by m − 2, the number of remaining nodes that can contribute edges in our accounting.The calculation establishes that A is strictly diagonally dominant and hence full rank.
Proof Let M ∈ {0, 1} ( m 2 )×m be the incidence matrix of a complete graph with m vertices.That is M has entry m e,v = 1 if vertex v occurs in edge e and 0 otherwise.Each row of M has two entries equal to 1; each column of M has m − 1 entries equal to 1.It is easy to see In general, it is easy to check that the matrix m×m matrix aI +b11 t has the eigenvector 1 with eigenvalue a + mb and m − 1 orthogonal eigenvectors with eigenvalue a.Note that each u i is perpendicular to 1. None of these eigenvectors is normalized to have length 1.Although the eigenvectors u i are certainly convenient, they are not unique.
To recover the eigenvectors of T t T , and hence those D t D, we can leverage the eigenvectors of M t M , which we know.The following generic observations are pertinent.If a matrix A has full SVD U SV t , then its transpose has full SVD A t = V SU t .As mentioned AA t and A t A share the same nontrivial eigenvalues.These can be recovered as the nontrivial diagonal entries of S 2 .Suppose we know the eigenvectors U of AA t = U S 2 U t .Since A t U = V S, then presumably we can recover some of the eigenvectors V as A t U S + , where S + is the diagonal pseudo-inverse of S. Applying the matrix inverse from before yields an explicit formula (with a and b defined as before): 1 ; z n = σ + ρD t P(Dx n ), a = 1 + ρp(c 2 + 1), b = 2ρc.

F.2 Steepest Descent
The updates x n+1 = x n − t n ∇h ρ (x n ) admit an exact solution for the line search parameter t n .Taking q n = ∇h ρ (x n ) as the gradient we have q n = (x n − u) + ρD t [Dx n − P ν (Dx n )], t n = q n 2 q n 2 + ρ Dq n 2 .

F.3 ADMM
Take y as the dual variable and λ as scaled multipliers.The formula for the MM algorithm applies in updating x n , except we replace ρ with µ and P(Dx n ) with y n − λ n : Multipliers follow the standard update.

F.4 Explicit Matrix Inverse
Both ADMM and MM reduce to solving a linear system.Fortunately, the Hessian for h ρ (x) reduces to a Householder-like matrix.First we note that it is trivial to multiply either C t or S t by a p 2 -vector.The more interesting problem is calculating (A     ( Bezanson et al., 2017).Additional packages used include Plots.jl(Breloff, 2021), GR.jl (Heinen et al., 2021), and (Udell et al., 2014).Numerical experiments were carried out on a Manjaro Linux 5.10.89-1desktop environment using 8 cores on an Intel 10900KF at 4.9 GHz and 32 GB RAM.

Figure 2 :
Figure 2: Sample images along the solution path of the search heuristic.Images are arranged from left to right as follows: reference, noisy input, and 90% reduction.

Figure 3 :
Figure 3: Example of a symmetric matrix X and its minimal representation x = trivec(X).

Table 2 :
Metric projection.Times are averaged over 3 replicates with standard deviations in parentheses.Reported iteration counts reflect the total inner iterations taken with outer iterations in parentheses.

Table 3 :
Convex regression.Times are reported as averages over 3 replicates.

Table 4 :
Convex clustering.Times reflect the total time spent generating candidate clusterings using Algorithm 2. Additional metrics correspond to the optimal clustering on the basis of maximal ARI.Time and clustering indices are averaged over 3 replicates with standard deviations reported in parentheses.

Table 5 :
Image denoising.Times reflect the total time spent generating candidate images, averaged over 3 replicates, ultimately achieving 90% reduction in total variation
t ρ A ρ ) −1 , where ∇h 2 ρ = A t ρ A ρ = I p + ρD t D.The reader can check the identitiesC t C = c 2 pI p and C t S = −c1 p 1 t p , S t C = −c1 p 1 tp and S t S = pI p .

Table 7 :
Performance of MM on convex regression using CG and LSQR.Both inner and outer iterations are reported with the latter in parantheses.It follows that A t ρ A ρ = (1 + ρp(c 2 + 1))I p − 2cρ1 p 1 t p .Applying the Sherman-Morrison formula to results in

Table 8 :
Performance of ADMM on convex regression using CG and LSQR.Both inner and outer iterations are reported with the latter in parantheses.